高斯勒让德(Gauss | 您所在的位置:网站首页 › python 二维积分 › 高斯勒让德(Gauss |
第四十四篇 高斯勒让德求解多重积分
多重积分
在工程分析中,经常需要在一个面积或体积上对函数进行积分。多重积分的解析方法在有限的情况下是可能的,但在这一篇中使用数值积分去求解。一维的函数积分详见重复牛顿-科特斯积分,重复高斯勒让德积分 考虑两个变量的函数在二维区域R上的积分,如下图所示。函数f(x, y)可以被认为是在R区域上以直角从平面上伸出的第三维空间。 在下图中,考虑函数f(x, y)在所示的矩形区域的积分。由于矩形的边界与笛卡尔坐标方向平行,变量之间没有关系,可以直接应用前面描述的积分方法。 例如,在每个方向上应用梯形规则将导致4个样本点(n = 4),如上图所示,矩形的每个角上都有一个样本点,从而产生该规则 在每个方向使用梯形法则去计算下面积分 在每个方向使用辛普森法则去计算下面的积分 高斯-勒让德法则也可以应用于这种类型的多重积分,但必须注意找到样本点的正确位置。一种方法是进行坐标变换,使每个方向上的积分极限变为±1。就可以使用这篇中提到得权值和采样点值。详细的表格对应数值可见高斯勒让德求积分 矩形区域的另一种方法是将样本点在每个坐标方向上的范围中点左右进行加权和比例系数相乘。 考虑下图所示的矩形区域上的两点高斯-勒让德积分。 在两个方向使用高斯-勒让德计算下面的积分 在每个方向使用3点高斯-勒让德去计算下面的积分 其中有一个主程序,和两个子程序,分别为高斯勒让德的样本点值和权值的子程序gauss_legendre;和局部坐标的形状函数的子程序fun_der。 #高斯拉盖尔法则 import numpy as np import B import math ndim=2;nsp=9 if ndim==1: nod=2 elif ndim==2: nod=4 elif ndim==3: nod=8 coord=np.zeros((nod,ndim));der=np.zeros((ndim,nod));fun=np.zeros(nod) coord[:,0]=(0.0,1.0,4.0,6.0);coord[:,1]=(0.0,2.0,3.0,1.0) samp=np.zeros((nsp,ndim));wt=np.zeros((nsp)) B.gauss_legendre(samp,wt) res=0 def f65(point): x=point[0];y=point[1]#;z=point[2] f65=x**2*y**2 return f65 for i in range(1,nsp+1): B.fun_der(fun,der,samp,i) res=res+np.linalg.det(np.dot(der,coord))*wt[i-1]*f65(np.dot(np.transpose(coord),fun)) print('高斯-勒让德法则的多重积分') print('维数',ndim) for i in range(1,nod+1): print('坐标 (x,y[,z])',coord[i-1,:]) print('样本点的数量',nsp) print('计算结果 ','{:13.4e}'.format(res)) gauss-legendre def gauss_legendre(samp,wt): nsp=samp.shape[0] ndim=samp.shape[1] if ndim==1: if nsp==1: samp[0,0]= 0.0 wt[0] = 2.0 elif nsp==2: samp[0,0]= -0.57735026918962576449 samp[1,0]= 0.57735026918962576449 wt[0]= 1.0 wt[1]= 1.0 elif nsp==3: samp[0,0]= -0.77459666924148337704 samp[1,0]= 0.0 samp[2,0]= 0.77459666924148337704 wt[0]= 0.55555555555555555556 wt[1]= 0.88888888888888888889 wt[2]= 0.55555555555555555556 elif nsp==4: samp[0,0]= -0.86113631159405257524 samp[1,0]= -0.33998104358485626481 samp[2,0]= 0.33998104358485626481 samp[3,0]= 0.86113631159405257524 wt[0]= 0.34785484513745385737 wt[1]= 0.65214515486254614271 wt[2]= 0.65214515486254614271 wt[3]= 0.34785484513745385737 elif nsp==5: samp[0,0]= -0.90617984593866399282 samp[1,0]= -0.53846931010568309105 samp[2,0]= 0.0 samp[3,0]= 0.53846931010568309105 samp[4,0]= 0.90617984593866399282 wt[0]= 0.23692688505618908751 wt[1]= 0.47862867049936646804 wt[2]= 0.56888888888888888889 wt[3]= 0.47862867049936646804 wt[4]= 0.23692688505618908751 else: print('积分点数量错误') elif ndim==2: if nsp==1: samp[0,0]= 0.0 samp[0,1]= 0.0 wt[0]= 4.0 elif nsp==4: samp[0,0]= -0.57735026918962576449 samp[1,0]= 0.57735026918962576449 samp[2,0]= -0.57735026918962576449 samp[3,0 |
今日新闻 |
推荐新闻 |
专题文章 |
CopyRight 2018-2019 实验室设备网 版权所有 |